Within Guy11

Differential peptide abundance

Perform differential abundance estimation using bootstrap and rank product methods

library(pepdiff)
d <- import_data(here::here("cleaned", "cleaned_data.csv"), 
                 gene_id = "molecule_list_name", 
                 treatment = "genotype", 
                 peptide = "peptide_modified_sequence", 
                 seconds = "timepoint")
## Rows: 20520 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): molecule_list_name, peptide_modified_sequence, genotype
## dbl (4): timepoint, bio_rep, tech_rep, total_area
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
comparisons <- data.frame(
  control = rep("Guy11", 5),
  c_seconds = rep(0, 5),
  treatment = rep("Guy11", 5),
  t_seconds = c(1,1.5,2,4,6)
)

many <- compare_many(d, comparisons, metrics = c("bootstrap_t", "rank_product"), iters=10000)
## Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done
many
## $`Guy11_1-Guy11_0`
## # A tibble: 285 × 31
##    gene_id     peptide     Guy11_1_1 Guy11_1_2 Guy11_1_3 Guy11_1.5_1 Guy11_1.5_2
##    <chr>       <chr>           <dbl>     <dbl>     <dbl>       <dbl>       <dbl>
##  1 gi|1456017… KLS[+80]AS… 0.000280  0.000330  0.000298    0.000745    0.000225 
##  2 gi|1456017… KLSASGS[+8… 0.00269   0.00348   0.00257     0.00285     0.00312  
##  3 gi|1456020… QTRT[+80]P… 0.00161   0.00121   0.00148     0.00243     0.000665 
##  4 gi|1456020… TPS[+80]PA… 0.00347   0.00198   0.00406     0.00356     0.00334  
##  5 gi|1456020… NAFLFGEVTT… 0.000497  0.000686  0.00111     0.000730    0.000967 
##  6 gi|1456023… AAEEDVTTT[… 0.0000318 0.0000318 0.0000318   0.0000318   0.0000318
##  7 gi|1456023… VQGNAYGGTG… 0.0000971 0.0000971 0.000290    0.000382    0.000526 
##  8 gi|1456026… AVGS[+80]R… 0.00362   0.00695   0.00838     0.00618     0.00766  
##  9 gi|1456026… AVGSRS[+80… 0.00370   0.00701   0.00834     0.00625     0.00761  
## 10 gi|1456026… DRS[+80]SG… 0.00195   0.000336  0.00170     0.00354     0.00107  
## # … with 275 more rows, and 24 more variables: Guy11_1.5_3 <dbl>,
## #   Guy11_0_1 <dbl>, Guy11_0_2 <dbl>, Guy11_0_3 <dbl>, fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_Guy11_1_1 <dbl>, unreplaced_Guy11_1_2 <dbl>,
## #   unreplaced_Guy11_1_3 <dbl>, unreplaced_Guy11_1.5_1 <dbl>,
## #   unreplaced_Guy11_1.5_2 <dbl>, unreplaced_Guy11_1.5_3 <dbl>,
## #   unreplaced_Guy11_0_1 <dbl>, unreplaced_Guy11_0_2 <dbl>, …
## 
## $`Guy11_1.5-Guy11_0`
## # A tibble: 285 × 25
##    gene_id    peptide    Guy11_1.5_1 Guy11_1.5_2 Guy11_1.5_3 Guy11_0_1 Guy11_0_2
##    <chr>      <chr>            <dbl>       <dbl>       <dbl>     <dbl>     <dbl>
##  1 gi|145601… KLS[+80]A…   0.000745    0.000225    0.000117  0.00108   0.000404 
##  2 gi|145601… KLSASGS[+…   0.00285     0.00312     0.00168   0.00215   0.00260  
##  3 gi|145602… QTRT[+80]…   0.00243     0.000665    0.000749  0.00287   0.00222  
##  4 gi|145602… TPS[+80]P…   0.00356     0.00334     0.00297   0.00128   0.00158  
##  5 gi|145602… NAFLFGEVT…   0.000730    0.000967    0.000520  0.000608  0.000626 
##  6 gi|145602… AAEEDVTTT…   0.0000318   0.0000318   0.0000318 0.0000318 0.0000318
##  7 gi|145602… VQGNAYGGT…   0.000382    0.000526    0.000272  0.0000971 0.0000971
##  8 gi|145602… AVGS[+80]…   0.00618     0.00766     0.00518   0.00193   0.00167  
##  9 gi|145602… AVGSRS[+8…   0.00625     0.00761     0.00518   0.00197   0.00167  
## 10 gi|145602… DRS[+80]S…   0.00354     0.00107     0.00190   0.000275  0.000275 
## # … with 275 more rows, and 18 more variables: Guy11_0_3 <dbl>,
## #   fold_change <dbl>, treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_Guy11_1.5_1 <dbl>, unreplaced_Guy11_1.5_2 <dbl>,
## #   unreplaced_Guy11_1.5_3 <dbl>, unreplaced_Guy11_0_1 <dbl>,
## #   unreplaced_Guy11_0_2 <dbl>, unreplaced_Guy11_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …
## 
## $`Guy11_2-Guy11_0`
## # A tibble: 285 × 25
##    gene_id  peptide  Guy11_2_1 Guy11_2_2 Guy11_2_3 Guy11_0_1 Guy11_0_2 Guy11_0_3
##    <chr>    <chr>        <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 gi|1456… KLS[+80… 0.000143   0.000174 0.000578  0.00108   0.000404  0.0000917
##  2 gi|1456… KLSASGS… 0.00155    0.00261  0.00122   0.00215   0.00260   0.00164  
##  3 gi|1456… QTRT[+8… 0.00233    0.00110  0.00135   0.00287   0.00222   0.000471 
##  4 gi|1456… TPS[+80… 0.00339    0.00377  0.00356   0.00128   0.00158   0.000209 
##  5 gi|1456… NAFLFGE… 0.000840   0.00107  0.000352  0.000608  0.000626  0.000501 
##  6 gi|1456… AAEEDVT… 0.0000318  0.000252 0.0000318 0.0000318 0.0000318 0.0000318
##  7 gi|1456… VQGNAYG… 0.0000971  0.00135  0.00111   0.0000971 0.0000971 0.0000971
##  8 gi|1456… AVGS[+8… 0.000938   0.00144  0.00358   0.00193   0.00167   0.00936  
##  9 gi|1456… AVGSRS[… 0.000938   0.00144  0.00358   0.00197   0.00167   0.00934  
## 10 gi|1456… DRS[+80… 0.00140    0.00130  0.00445   0.000275  0.000275  0.000275 
## # … with 275 more rows, and 17 more variables: fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_Guy11_2_1 <dbl>, unreplaced_Guy11_2_2 <dbl>,
## #   unreplaced_Guy11_2_3 <dbl>, unreplaced_Guy11_0_1 <dbl>,
## #   unreplaced_Guy11_0_2 <dbl>, unreplaced_Guy11_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …
## 
## $`Guy11_4-Guy11_0`
## # A tibble: 285 × 25
##    gene_id  peptide  Guy11_4_1 Guy11_4_2 Guy11_4_3 Guy11_0_1 Guy11_0_2 Guy11_0_3
##    <chr>    <chr>        <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 gi|1456… KLS[+80… 0.00125    0.000264  0.000475 0.00108   0.000404  0.0000917
##  2 gi|1456… KLSASGS… 0.00298    0.00202   0.00228  0.00215   0.00260   0.00164  
##  3 gi|1456… QTRT[+8… 0.00184    0.00116   0.000869 0.00287   0.00222   0.000471 
##  4 gi|1456… TPS[+80… 0.00514    0.00371   0.00457  0.00128   0.00158   0.000209 
##  5 gi|1456… NAFLFGE… 0.00184    0.00153   0.00166  0.000608  0.000626  0.000501 
##  6 gi|1456… AAEEDVT… 0.0000318  0.000194  0.000111 0.0000318 0.0000318 0.0000318
##  7 gi|1456… VQGNAYG… 0.00449    0.00400   0.00538  0.0000971 0.0000971 0.0000971
##  8 gi|1456… AVGS[+8… 0.0196     0.0128    0.00810  0.00193   0.00167   0.00936  
##  9 gi|1456… AVGSRS[… 0.0197     0.0128    0.00818  0.00197   0.00167   0.00934  
## 10 gi|1456… DRS[+80… 0.00845    0.00328   0.00431  0.000275  0.000275  0.000275 
## # … with 275 more rows, and 17 more variables: fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_Guy11_4_1 <dbl>, unreplaced_Guy11_4_2 <dbl>,
## #   unreplaced_Guy11_4_3 <dbl>, unreplaced_Guy11_0_1 <dbl>,
## #   unreplaced_Guy11_0_2 <dbl>, unreplaced_Guy11_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …
## 
## $`Guy11_6-Guy11_0`
## # A tibble: 285 × 25
##    gene_id  peptide  Guy11_6_1 Guy11_6_2 Guy11_6_3 Guy11_0_1 Guy11_0_2 Guy11_0_3
##    <chr>    <chr>        <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 gi|1456… KLS[+80… 0.0000917  0.000129  0.000127 0.00108   0.000404  0.0000917
##  2 gi|1456… KLSASGS… 0.000851   0.00129   0.000982 0.00215   0.00260   0.00164  
##  3 gi|1456… QTRT[+8… 0.00264    0.000487  0.00362  0.00287   0.00222   0.000471 
##  4 gi|1456… TPS[+80… 0.00430    0.00227   0.00558  0.00128   0.00158   0.000209 
##  5 gi|1456… NAFLFGE… 0.000588   0.00185   0.000952 0.000608  0.000626  0.000501 
##  6 gi|1456… AAEEDVT… 0.000450   0.000128  0.000647 0.0000318 0.0000318 0.0000318
##  7 gi|1456… VQGNAYG… 0.00184    0.00203   0.00158  0.0000971 0.0000971 0.0000971
##  8 gi|1456… AVGS[+8… 0.00123    0.0111    0.000650 0.00193   0.00167   0.00936  
##  9 gi|1456… AVGSRS[… 0.00123    0.0111    0.000650 0.00197   0.00167   0.00934  
## 10 gi|1456… DRS[+80… 0.00335    0.00211   0.00349  0.000275  0.000275  0.000275 
## # … with 275 more rows, and 17 more variables: fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_Guy11_6_1 <dbl>, unreplaced_Guy11_6_2 <dbl>,
## #   unreplaced_Guy11_6_3 <dbl>, unreplaced_Guy11_0_1 <dbl>,
## #   unreplaced_Guy11_0_2 <dbl>, unreplaced_Guy11_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …

Plot result heatmaps

lgd <- ComplexHeatmap::Legend(direction = "horizontal",
  col_fun = circlize::colorRamp2(
                           seq(-2, 2, length.out = 11),
                           rev(RColorBrewer::brewer.pal(11, "RdBu"))
                         ),
  title = "Log 2 Fold Change")


plot_heatmap(many, metric = "bootstrap_t", log = TRUE)#,  col_order = col_order)
## Joining, by = "comparison"
## Joining, by = "gene_peptide"
ComplexHeatmap::draw(lgd, x = grid::unit(8, "in"), y = grid::unit(1, "in"))

plot_heatmap(many, metric = "rank_product", log = TRUE)
## Joining, by = "comparison"
## Joining, by = "gene_peptide"
ComplexHeatmap::draw(lgd, x = grid::unit(8, "in"), y = grid::unit(1, "in"))

#add colour even when the point isn't significant
plot_heatmap(many, metric = "bootstrap_t", log = TRUE)#,  col_order = col_order, only_sig_points = FALSE)
## Joining, by = "comparison"
## Joining, by = "gene_peptide"
ComplexHeatmap::draw(lgd, x = grid::unit(8, "in"), y = grid::unit(1, "in"))

volcano_plot(many, metric = "bootstrap_t")

writexl::write_xlsx(many, here::here("analysis","003_differential_abundance_Guy11_w_geno_estimates.xlsx"))

Within dpmk1

comparisons <- data.frame(
  control = rep("dpmk1", 5),
  c_seconds = rep(0, 5),
  treatment = rep("dpmk1", 5),
  t_seconds = c(1,1.5,2,4,6)
)

many <- compare_many(d, comparisons, metrics = c("bootstrap_t", "rank_product"), iters=10000)
## Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done  Rank Product analysis for unpaired case 
##  
## 
##  done
many
## $`dpmk1_1-dpmk1_0`
## # A tibble: 285 × 31
##    gene_id     peptide     dpmk1_1_1 dpmk1_1_2 dpmk1_1_3 dpmk1_1.5_1 dpmk1_1.5_2
##    <chr>       <chr>           <dbl>     <dbl>     <dbl>       <dbl>       <dbl>
##  1 gi|1456017… KLS[+80]AS… 0.000308  0.000281  0.000335    0.000456    0.000157 
##  2 gi|1456017… KLSASGS[+8… 0.00202   0.00157   0.00455     0.00226     0.00200  
##  3 gi|1456020… QTRT[+80]P… 0.000876  0.00201   0.00139     0.00113     0.000698 
##  4 gi|1456020… TPS[+80]PA… 0.00240   0.00240   0.00231     0.00334     0.00178  
##  5 gi|1456020… NAFLFGEVTT… 0.00102   0.00123   0.00317     0.000900    0.000719 
##  6 gi|1456023… AAEEDVTTT[… 0.0000318 0.0000318 0.0000318   0.0000318   0.0000318
##  7 gi|1456023… VQGNAYGGTG… 0.000107  0.000365  0.000286    0.000351    0.000645 
##  8 gi|1456026… AVGS[+80]R… 0.00422   0.00564   0.0171      0.0110      0.00557  
##  9 gi|1456026… AVGSRS[+80… 0.00424   0.00566   0.0172      0.0110      0.00561  
## 10 gi|1456026… DRS[+80]SG… 0.000721  0.00150   0.000447    0.00155     0.00257  
## # … with 275 more rows, and 24 more variables: dpmk1_1.5_3 <dbl>,
## #   dpmk1_0_1 <dbl>, dpmk1_0_2 <dbl>, dpmk1_0_3 <dbl>, fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_dpmk1_1_1 <dbl>, unreplaced_dpmk1_1_2 <dbl>,
## #   unreplaced_dpmk1_1_3 <dbl>, unreplaced_dpmk1_1.5_1 <dbl>,
## #   unreplaced_dpmk1_1.5_2 <dbl>, unreplaced_dpmk1_1.5_3 <dbl>,
## #   unreplaced_dpmk1_0_1 <dbl>, unreplaced_dpmk1_0_2 <dbl>, …
## 
## $`dpmk1_1.5-dpmk1_0`
## # A tibble: 285 × 25
##    gene_id    peptide    dpmk1_1.5_1 dpmk1_1.5_2 dpmk1_1.5_3 dpmk1_0_1 dpmk1_0_2
##    <chr>      <chr>            <dbl>       <dbl>       <dbl>     <dbl>     <dbl>
##  1 gi|145601… KLS[+80]A…   0.000456    0.000157    0.0000917 0.00123   0.000316 
##  2 gi|145601… KLSASGS[+…   0.00226     0.00200     0.00160   0.000732  0.000321 
##  3 gi|145602… QTRT[+80]…   0.00113     0.000698    0.000195  0.00304   0.00271  
##  4 gi|145602… TPS[+80]P…   0.00334     0.00178     0.000529  0.000498  0.000462 
##  5 gi|145602… NAFLFGEVT…   0.000900    0.000719    0.000654  0.000306  0.000484 
##  6 gi|145602… AAEEDVTTT…   0.0000318   0.0000318   0.0000318 0.0000318 0.0000318
##  7 gi|145602… VQGNAYGGT…   0.000351    0.000645    0.000387  0.0000971 0.0000971
##  8 gi|145602… AVGS[+80]…   0.0110      0.00557     0.0102    0.0297    0.0106   
##  9 gi|145602… AVGSRS[+8…   0.0110      0.00561     0.0102    0.0296    0.0105   
## 10 gi|145602… DRS[+80]S…   0.00155     0.00257     0.000440  0.000275  0.000275 
## # … with 275 more rows, and 18 more variables: dpmk1_0_3 <dbl>,
## #   fold_change <dbl>, treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_dpmk1_1.5_1 <dbl>, unreplaced_dpmk1_1.5_2 <dbl>,
## #   unreplaced_dpmk1_1.5_3 <dbl>, unreplaced_dpmk1_0_1 <dbl>,
## #   unreplaced_dpmk1_0_2 <dbl>, unreplaced_dpmk1_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …
## 
## $`dpmk1_2-dpmk1_0`
## # A tibble: 285 × 25
##    gene_id  peptide  dpmk1_2_1 dpmk1_2_2 dpmk1_2_3 dpmk1_0_1 dpmk1_0_2 dpmk1_0_3
##    <chr>    <chr>        <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 gi|1456… KLS[+80… 0.000251  0.0000917 0.000352  0.00123   0.000316  0.0000917
##  2 gi|1456… KLSASGS… 0.00145   0.00197   0.00464   0.000732  0.000321  0.000321 
##  3 gi|1456… QTRT[+8… 0.00123   0.000270  0.00171   0.00304   0.00271   0.000242 
##  4 gi|1456… TPS[+80… 0.00246   0.00262   0.00375   0.000498  0.000462  0.00138  
##  5 gi|1456… NAFLFGE… 0.000849  0.000837  0.00234   0.000306  0.000484  0.000206 
##  6 gi|1456… AAEEDVT… 0.0000318 0.0000318 0.0000318 0.0000318 0.0000318 0.0000318
##  7 gi|1456… VQGNAYG… 0.000275  0.00370   0.00145   0.0000971 0.0000971 0.0000971
##  8 gi|1456… AVGS[+8… 0.00872   0.00601   0.0144    0.0297    0.0106    0.00329  
##  9 gi|1456… AVGSRS[… 0.00866   0.00606   0.0144    0.0296    0.0105    0.00329  
## 10 gi|1456… DRS[+80… 0.00106   0.00491   0.00193   0.000275  0.000275  0.000275 
## # … with 275 more rows, and 17 more variables: fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_dpmk1_2_1 <dbl>, unreplaced_dpmk1_2_2 <dbl>,
## #   unreplaced_dpmk1_2_3 <dbl>, unreplaced_dpmk1_0_1 <dbl>,
## #   unreplaced_dpmk1_0_2 <dbl>, unreplaced_dpmk1_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …
## 
## $`dpmk1_4-dpmk1_0`
## # A tibble: 285 × 25
##    gene_id  peptide  dpmk1_4_1 dpmk1_4_2 dpmk1_4_3 dpmk1_0_1 dpmk1_0_2 dpmk1_0_3
##    <chr>    <chr>        <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 gi|1456… KLS[+80…  0.000447 0.000671   0.000263 0.00123   0.000316  0.0000917
##  2 gi|1456… KLSASGS…  0.00276  0.00395    0.00294  0.000732  0.000321  0.000321 
##  3 gi|1456… QTRT[+8…  0.000959 0.000550   0.000818 0.00304   0.00271   0.000242 
##  4 gi|1456… TPS[+80…  0.00407  0.00296    0.00310  0.000498  0.000462  0.00138  
##  5 gi|1456… NAFLFGE…  0.00306  0.00160    0.00281  0.000306  0.000484  0.000206 
##  6 gi|1456… AAEEDVT…  0.000367 0.0000750  0.000265 0.0000318 0.0000318 0.0000318
##  7 gi|1456… VQGNAYG…  0.00607  0.00776    0.00710  0.0000971 0.0000971 0.0000971
##  8 gi|1456… AVGS[+8…  0.0180   0.0160     0.0178   0.0297    0.0106    0.00329  
##  9 gi|1456… AVGSRS[…  0.0179   0.0161     0.0178   0.0296    0.0105    0.00329  
## 10 gi|1456… DRS[+80…  0.00375  0.00755    0.00298  0.000275  0.000275  0.000275 
## # … with 275 more rows, and 17 more variables: fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_dpmk1_4_1 <dbl>, unreplaced_dpmk1_4_2 <dbl>,
## #   unreplaced_dpmk1_4_3 <dbl>, unreplaced_dpmk1_0_1 <dbl>,
## #   unreplaced_dpmk1_0_2 <dbl>, unreplaced_dpmk1_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …
## 
## $`dpmk1_6-dpmk1_0`
## # A tibble: 285 × 25
##    gene_id  peptide  dpmk1_6_1 dpmk1_6_2 dpmk1_6_3 dpmk1_0_1 dpmk1_0_2 dpmk1_0_3
##    <chr>    <chr>        <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 gi|1456… KLS[+80… 0.000101  0.0000992 0.000278  0.00123   0.000316  0.0000917
##  2 gi|1456… KLSASGS… 0.00163   0.000856  0.00132   0.000732  0.000321  0.000321 
##  3 gi|1456… QTRT[+8… 0.000110  0.00323   0.000200  0.00304   0.00271   0.000242 
##  4 gi|1456… TPS[+80… 0.00130   0.00591   0.00147   0.000498  0.000462  0.00138  
##  5 gi|1456… NAFLFGE… 0.00108   0.000803  0.00335   0.000306  0.000484  0.000206 
##  6 gi|1456… AAEEDVT… 0.0000318 0.000391  0.0000318 0.0000318 0.0000318 0.0000318
##  7 gi|1456… VQGNAYG… 0.00312   0.00126   0.00280   0.0000971 0.0000971 0.0000971
##  8 gi|1456… AVGS[+8… 0.0160    0.00160   0.0290    0.0297    0.0106    0.00329  
##  9 gi|1456… AVGSRS[… 0.0158    0.00162   0.0290    0.0296    0.0105    0.00329  
## 10 gi|1456… DRS[+80… 0.00305   0.00360   0.00451   0.000275  0.000275  0.000275 
## # … with 275 more rows, and 17 more variables: fold_change <dbl>,
## #   treatment_mean_count <dbl>, control_mean_count <dbl>,
## #   unreplaced_dpmk1_6_1 <dbl>, unreplaced_dpmk1_6_2 <dbl>,
## #   unreplaced_dpmk1_6_3 <dbl>, unreplaced_dpmk1_0_1 <dbl>,
## #   unreplaced_dpmk1_0_2 <dbl>, unreplaced_dpmk1_0_3 <dbl>,
## #   treatment_replicates <int>, control_replicates <int>,
## #   bootstrap_t_p_val <dbl>, bootstrap_t_fdr <dbl>, rank_prod_p2_p_val <dbl>, …

Plot result heatmaps

lgd <- ComplexHeatmap::Legend(direction = "horizontal",
  col_fun = circlize::colorRamp2(
                           seq(-2, 2, length.out = 11),
                           rev(RColorBrewer::brewer.pal(11, "RdBu"))
                         ),
  title = "Log 2 Fold Change")


plot_heatmap(many, metric = "bootstrap_t", log = TRUE)#,  col_order = col_order)
## Joining, by = "comparison"
## Joining, by = "gene_peptide"
ComplexHeatmap::draw(lgd, x = grid::unit(8, "in"), y = grid::unit(1, "in"))

plot_heatmap(many, metric = "rank_product", log = TRUE)
## Joining, by = "comparison"
## Joining, by = "gene_peptide"
ComplexHeatmap::draw(lgd, x = grid::unit(8, "in"), y = grid::unit(1, "in"))

#add colour even when the point isn't significant
plot_heatmap(many, metric = "bootstrap_t", log = TRUE)#,  col_order = col_order, only_sig_points = FALSE)
## Joining, by = "comparison"
## Joining, by = "gene_peptide"
ComplexHeatmap::draw(lgd, x = grid::unit(8, "in"), y = grid::unit(1, "in"))

volcano_plot(many, metric = "bootstrap_t")

writexl::write_xlsx(many, here::here("analysis","003_differential_abundance_dpmk_w_geno_estimates.xlsx"))